home *** CD-ROM | disk | FTP | other *** search
/ Business Shareware / Business Shareware.iso / start / hobby / aa_51 / osaturn.c < prev    next >
C/C++ Source or Header  |  1993-02-13  |  7KB  |  294 lines
  1. /* Elements of the orbit of the planet Saturn
  2.  * using formulas given by Meeus.
  3.  *
  4.  * This program seems to have much larger errors than 
  5.  * the ones for the other planets.  Error of 90 arc seconds
  6.  * in R.A. and 15 arc seconds in Dec observed re 1986 tables.
  7.  */
  8.  
  9. #include "planet.h"
  10. #include "kep.h"
  11. static double sinpsi = 0.0;
  12. static double cospsi = 0.0;
  13. static double sin2psi = 0.0;
  14. static double cos2psi = 0.0;
  15. static double sin3psi = 0.0;
  16. static double cos3psi = 0.0;
  17.  
  18. int pjup(), esat(), asat();
  19.  
  20. int osaturn(e,J)
  21. struct orbit *e;
  22. double J;
  23. {
  24. double p, q;
  25.  
  26. e->epoch = J;
  27. manoms(J);
  28. pjup();
  29. sinpsi = sin(psi);
  30. cospsi = cos(psi);
  31. sin2psi = 2.0*sinpsi*cospsi;
  32. cos2psi = cospsi*cospsi - sinpsi*sinpsi;
  33. sin3psi = sinpsi*cos2psi + cospsi*sin2psi;
  34. cos3psi = cospsi*cos2psi - sinpsi*sin2psi;
  35. e->M = M6;
  36.  
  37. #if OFDATE
  38. /* expansions for equinox of date */
  39. e->equinox = J;
  40. /* inclination */
  41. e->i = (( 0.00000004*T - 0.00001549)*T - 0.0039189)*T + 2.492519;
  42. /* arg perihelion */
  43. f = (( 0.00000992*T + 0.00097854)*T + 1.0852207)*T + 338.307800;
  44. e->w = mod360(f);
  45. /* ascending node */
  46. f = (( -0.00000531*T - 0.00015218)*T + 0.8731951)*T + 112.790414;
  47. e->W = mod360(f);
  48. /* mean longitude */
  49. f = (( -0.0000058*T + 0.0003245)*T + 1223.509884)*T + 266.564377;
  50. e->L = mod360(f);
  51. #else
  52. /* J2000 coordinates */
  53. e->equinox = J2000;
  54. /* inclination */
  55. e->i = (( 2e-9*T - 5.017e-5)*T + 0.0024449)*T + 2.486204;
  56. /* arg perihelion */
  57. f = (( 6.177e-6*T + 0.00070747)*T + 0.8220515)*T + 338.571353;
  58. e->w = mod360(f);
  59. /* ascending node */
  60. f = (( -1.589e-6*T - 0.00018997)*T - 0.2599254)*T + 113.923406;
  61. e->W = mod360(f);
  62. /* mean longitude */
  63. f = e->M + e->w + e->W;
  64. e->L = mod360(f);
  65. #endif
  66.  
  67. /* eccentricity */
  68. e->ecc = (( 0.00000000074*T - 0.000000728)*T - 0.00034550)*T + 0.05589232;
  69. /* mean distance */
  70. e->a = 9.554747;
  71.         
  72. /* perturbations in the mean longitude */
  73. p = ((0.016714*nu + 0.018150)*nu - 0.814181)*sinV
  74.     + (( -0.004100*nu + 0.160906)*nu - 0.010497)*cosV
  75.     + 0.007581*sin2V
  76.     - 0.007986*sinW;
  77. p = p - 0.148811*sinze
  78.      - 0.040786*sin2ze
  79.      - 0.015208*sin3ze
  80.      - 0.006339*sin4ze;
  81. f = -0.006244
  82.     + (0.008931 + 0.002728*nu)*sinze
  83.     - 0.016500*sin2ze
  84.     - 0.005775*sin3ze
  85.     + (0.081344 + 0.003206*nu)*cosze
  86.     + 0.015019*cos2ze;
  87. p += f*sinQ;
  88.  
  89. f =  (0.085581 + 0.002494*nu)*sinze
  90.     + (0.025328 - 0.003117*nu)*cosze
  91.     + 0.014394*cos2ze
  92.     + 0.006319*cos3ze;
  93. p += f*cosQ;
  94.  
  95. f =  0.006369*sinze
  96.     + 0.009156*sin2ze
  97.     + 0.007525*sin3psi;
  98. p += f*sin2Q;
  99.  
  100. f = -0.005236*cosze
  101.     - 0.007736*cos2ze
  102.     - 0.007528*cos3psi;
  103. p += f*cos2Q;        
  104. e->L += p;
  105.  
  106. /* perturbations in the perihelion */
  107. q = ((-0.001533*nu + 0.007186)*nu + 0.077108)*sinV
  108.    + ((-0.000536*nu - 0.014766)*nu + 0.045803)*cosV
  109.     - 0.007075*sinze;
  110.  
  111. f = -0.075825*sinze
  112.     - 0.024839*sin2ze
  113.     - 0.008631*sin3ze;
  114. q += f*sinQ;
  115.  
  116. f = -0.072586
  117.     - 0.150383*cosze
  118.     + 0.026897*cos2ze
  119.     + 0.010053*cos3ze;
  120. q += f*cosQ;
  121.  
  122. f = -(0.013597 +0.001719*nu)*sinze
  123.     + (-0.007742 + 0.001517*nu)*cosze
  124.     + (0.013586 - 0.001375*nu)*cos2ze;
  125. q += f*sin2Q;
  126.  
  127. f =  (-0.013667 + 0.001239*nu)*sinze
  128.     + 0.011981*sin2ze
  129.     + (0.014861 + 0.001136*nu)*cosze
  130.     - (0.013064 + 0.001628*nu)*cos2ze;
  131. q += f*cos2Q;        
  132. e->w += q;
  133.  
  134. f = p - q/(e->ecc);
  135. e->M += f;
  136. esat(e);
  137. asat(e);
  138. return(0);
  139. }
  140.  
  141.  
  142. int esat(e)
  143. struct orbit *e;
  144. {
  145. double sin4Q, cos4Q;
  146. double q;
  147.  
  148. /* perturbations in the eccentricity */
  149. q =    ((0.0000091*nu + 0.0002548)*nu - 0.0007927)*sinV;
  150. q += (( -0.0000253*nu + 0.0001226)*nu + 0.0013381)*cosV;
  151. q +=  ( -0.0000121*nu + 0.0000248)*sin2V;
  152. q -=    (0.0000091*nu + 0.0000305)*cos2V;
  153. q += 0.0000412*sin2ze;
  154. f =  0.0012415
  155.     + ( 0.0000390 - 0.0000617*nu)*sinze
  156.     + ( 0.0000165 - 0.0000204*nu)*sin2ze;
  157. f = f + 0.0026599*cosze
  158.       - 0.0004687*cos2ze
  159.       - 0.0001870*cos3ze
  160.       - 0.0000821*cos4ze
  161.       - 0.0000377*cos5ze;
  162. f = f + 0.0000497*cos2psi;
  163. q += f*sinQ;
  164.  
  165. f =   0.0000163 - 0.0000611*nu
  166.     - 0.0012696*sinze
  167.     - 0.0004200*sin2ze
  168.     - 0.0001503*sin3ze
  169.     - 0.0000619*sin4ze
  170.     - 0.0000268*sin5ze;
  171. f = f - ( 0.0000282 + 0.0001306*nu)*cosze
  172.       + (-0.0000086 + 0.0000230*nu)*cos2ze;
  173. f = f   + 0.0000461*sin2psi;
  174. q += f*cosQ;
  175.  
  176. f =  -0.0000350
  177.     + (0.0002211 - 0.0000286*nu)*sinze
  178.     - 0.0002208*sin2ze
  179.     - 0.0000568*sin3ze
  180.     - 0.0000346*sin4ze;
  181. f = f - (0.0002780 + 0.0000222*nu)*cosze
  182.       + (0.0002022 + 0.0000263*nu)*cos2ze
  183.       + 0.0000248*cos3ze
  184.       + 0.0000242*sin3psi
  185.       + 0.0000467*cos3psi;
  186. q += f*sin2Q;
  187.  
  188. f =   -0.0000490
  189.     - ( 0.0002842 + 0.0000279*nu)*sinze
  190.     + ( 0.0000128 + 0.0000226*nu)*sin2ze;
  191. f = f + 0.0000224*sin3ze
  192.     + (-0.0001594 + 0.0000282*nu)*cosze
  193.     + ( 0.0002162 - 0.0000207*nu)*cos2ze;
  194. f = f + 0.0000561*cos3ze
  195.       + 0.0000343*cos4ze
  196.       + 0.0000469*sin3psi
  197.       - 0.0000242*cos3psi;
  198. q += f*cos2Q;
  199.  
  200. f =  -0.0000205*sinze
  201.     + 0.0000262*sin3ze;
  202. q += f*sin3Q;
  203.  
  204. f =   0.0000208*cosze
  205.     - 0.0000271*cos3ze;
  206. q += f*cos3Q;
  207.  
  208. sin4Q = 2.0*sin2Q*cos2Q;
  209. cos4Q = cos2Q*cos2Q - sin2Q*sin2Q;
  210. q = q - 0.0000382*cos3ze*sin4Q
  211.       - 0.0000376*sin3ze*cos4Q;
  212. e->ecc += q;
  213. return(0);
  214. }
  215.  
  216.  
  217.  
  218. int asat(e)
  219. struct orbit *e;
  220. {
  221. double q;
  222.  
  223. /* perturbations in the semimajor axis */
  224. q =   0.000572*nu*sinV
  225.     + 0.002933*cosV
  226.     + 0.033629*cosze
  227.     - 0.003081*cos2ze
  228.     - 0.001423*cos3ze
  229.     - 0.000671*cos4ze
  230.     - 0.000320*cos5ze;
  231.  
  232. f =   0.001098
  233.     - 0.002812*sinze
  234.     + 0.000688*sin2ze
  235.     - 0.000393*sin3ze
  236.     - 0.000228*sin4ze
  237.     + 0.002138*cosze
  238.     - 0.000999*cos2ze
  239.     - 0.000642*cos3ze
  240.     - 0.000325*cos4ze;
  241. q += f*sinQ;
  242.  
  243. f =  -0.000890
  244.     + 0.002206*sinze
  245.     - 0.001590*sin2ze
  246.     - 0.000647*sin3ze
  247.     - 0.000344*sin4ze
  248.     + 0.002885*cosze
  249.     + ( 0.002172 + 0.000102*nu)*cos2ze
  250.     + 0.000296*cos3ze;
  251. q += f*cosQ;
  252.  
  253. f =  -0.000267*sin2ze
  254.     - 0.000778*cosze
  255.     + 0.000495*cos2ze
  256.     + 0.000250*cos3ze;
  257. q += f*sin2Q;
  258.  
  259. f =  -0.000856*sinze
  260.     + 0.000441*sin2ze
  261.     + 0.000296*cos2ze
  262.     + 0.000211*cos3ze;
  263. q += f*cos2Q;
  264.  
  265. f =  -0.000427*sinze
  266.     + 0.000398*sin3ze;
  267. q += f*sin3Q;
  268. f =   0.000344*cosze
  269.     - 0.000427*cos3ze;
  270. q += f*cos3Q;       
  271. e->a += q;
  272. return(0);
  273. }
  274.  
  275.  
  276. /* Corrections to the position to be applied
  277.  * after solving Kepler's equation
  278.  */
  279. int csaturn(e)
  280. struct orbit *e;
  281. {
  282. double p;
  283.  
  284. /* perturbations to the heliocentric latitude */
  285. f = 0.000747*sinQ +0.001069*cosQ;
  286. p = f*cosze;
  287. f = 0.002108*sin2ze +0.001261*cos2ze;
  288. p = p + f*sin2Q;
  289. f = 0.001236*sin2ze -0.002075*cos2ze;
  290. p = p + f*cos2Q;
  291. e->plat = p*DTR;
  292. return(0);
  293. }
  294.